% REDCONF 
% Function testing the amount of red noise in a signal, and calculating the
% confidence levels using three methods : 
% -based on the spectral analysis of nsim simulations of red noise, using
% Monte-Carlo algorythme
% -based on a chi2 test performed on the mean spectrum of the nsim red noise
% calculated
% -based on the chi2 test performed on the theoretical estimate of the red
% noise spectrum.
% The spectral analysis method used for analysing the data and the red noise 
% is the MTM method, using Matlab function pmtm.
%
% The sample step must be constant, implying a re-sampling of the data
% prior to the analysis. Pre-processing of the data such as detrending, if applied, must
% also be performed prior to the analysis.
%
% Redconf function calls the subroutines rhoAR1.m, Rednoise.m and conflevel.m
% 
%
% Inputs : 
% timex= one column array with time or depth corresponding to each
% values of the signal studied
% datax= one column array with each values of the signal studied
% dt= sample step
% nsim= number of red noise simulations performed within the Monte Carlo loop.
% 1000 or 1500 are usually enough.
% nw= number of tapers used by the MTM method.

% Outputs :
% rho= lag-1 autocorrelation coefficient.
% Msepcred= 1 column array containing the mean spectrum of the nsim red noise modelings
% specred= n frequencies/nsim array containing all the spectra calculated
% for each red noise estimation.
% po= 1 column array containing the power values of the data spectrum
% fd= 1 column array containing the frequency values of the data spectrum
% pr= 1 column array containing the power values of the mean spectrum of the 
% read noise 
% fr= 1 column array containing the frequency values of the red noise mean 
% spectrum (fr should be equal to fd)
% 
% Doroth�e Husson, November 2012

function [Mspecred,specred,po,fd,fr,tabMC,tabchi,tabtchi,theored,rho]=redconf(datax,timex,dt,nsim,nw)

% calculation of the lag-1 autocorrelation coefficient
[rho]=rhoAR1(datax);

% calculation of nsim red noise models
[redtab,nt]=Rednoise(timex,rho,nsim);

% spectral analysis of the data using MTM method (matlab pmtm function)
% nw is the time-bandwidth product for the discrete prolate spheroidal sequences used as data tapers.
% If you specify nw as the empty vector [], a default value of 4 is used. 
% Other typical choices are 2, 5/2, 3, or 7/2 (from Matlab documentation)

% il vaut mieux enlever la valeur moyenne de tous les signaux
% (donnees reelles et bruits simules) avant de calculer les pmtm
% pour mieux comparer les aires
% normalization of the data;
datan=[];
i=1;
nd=length(datax);
dataxm=sum(datax)/nd;
datan=datax-dataxm;

% spectral analysis of the data
[po,w]=pmtm(datan,nw);
fd=w/(2*pi*dt);
npo=length(po);

% calculation of the area of the data power spectrum
Ax=(sum(po))/npo;

% spectral analysis of the nsim red noise signals and power normalization of
% each specra.
red2=[];
specred=[];
i=1;

for i=1:nsim
    red2=redtab(:,i); 
    red2m=sum(red2)/nt;
    red2n=red2-red2m; % cf commentaire precedent (enlever la moyenne de red2)
    [pr,wr]=pmtm(red2n,nw);
    npr=size(pr,1);
    %calculation of the area of the power spectrum of the rednoise
    Ar=sum(pr)/npr;
    pr=pr*(Ax/Ar);
    specred(:,i)=pr;
end
% Calculation of the mean rednoise spectra and normalization of the power.
i=1;
j=1;
Mpr=0;
Mspecred=[];
for i=1:npr
    Mpr=specred(i,1);
    for j=2:nsim
        Mpr=Mpr+specred(i,j);
    end
    Mspecred(i)=Mpr/nsim;
end
fr=wr/(2*pi*dt);

% calculation of the confidence levels
[tabMC,tabchi,tabtchi,theored]=conflevel(Mspecred,specred,nsim,nw,Ax,npr,dt,rho,fr);
end
